REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 
regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 

Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


5c.  PROGRAM  ELEMENT  NUMBER 

611102 


5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  Public  Release;  Distribution  Unlimited 

13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 

14.  ABSTRACT 

Freezing  and  thawing  of  soils  are  processes  that  span  multiple  spatial  scales  and  involve  multi-physics  phenomena 
such  as  phase  change  and  multiphase  flow.  Of  particular  importance  is  freezing  and  thawing  in  soils  susceptible  to 
heaving.  Frost  heave  and  thawing  of  ice-rich  soils  directly  affect  Army  operations  in  cold  regions  and  the  state  of 
infrastructure.  This  research  was  focused  on  the  study  of  the  fundamental  processes  and  causes  of  frost  heaving  and 
thaw-weakening  in  soils.  A  model  describing  frost  heaving  in  soils  has  been  developed  that  well  captures  the 

15.  SUBJECT  TERMS 

Freeze-thaw  cycle,  frost  heave,  ice  growth,  ice  lensing,  phase  change,  soil  thawing,  thaw-weakening,  mass  transfer,  heat  transfer, 
frost  susceptibility,  constitutive  model 

17.  LIMITATION  OF  15.  NUMBER 

ABSTRACT  OF  PAGES 

UU 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39. 1 8 


19a.  NAME  OF  RESPONSIBLE  PERSON 

Radoslaw  Michalowski _ 

19b.  TELEPHONE  NUMBER 
734-763-2146 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

UU 

UU 

UU 

7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

University  of  Michigan  -  Ann  Arbor 
Regents  of  the  University  of  Michigan 
3003  S.  State  St 

Ann  Arbor,  MI  48109  -1274 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND 
ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 

1 1 .  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

53315-EV.14 


5d.  PROJECT  NUMBER 


3.  DATES  COVERED  (From  -  To) 

1 -Feb-2009  -  31 -Jan-20 13 

5a.  CONTRACT  NUMBER 

W91  INF-08-1-0376 _ 

5b.  GRANT  NUMBER 


2.  REPORT  TYPE 

Final  Report 

4.  TITLE  AND  SUBTITLE 

Multi-Scale  Process  of  soil  Freezing,  Thawing,  and 
Thaw- Settlement 


6.  AUTHORS 

Radoslaw  L.  Michalowski  and  Yao  Zhang 


1.  REPORT  DATE  (DD-MM-YYYY) 
18-02-2013 


Report  Title 

Multi-Scale  Process  of  soil  Freezing,  Thawing,  and  Thaw-Settlement 

ABSTRACT 

Freezing  and  thawing  of  soils  are  processes  that  span  multiple  spatial  scales  and  involve  multi-physics  phenomena  such  as  phase  change  and 
multiphase  flow.  Of  particular  importance  is  freezing  and  thawing  in  soils  susceptible  to  heaving.  Frost  heave  and  thawing  of  ice-rich  soils 
directly  affect  Army  operations  in  cold  regions  and  the  state  of  infrastructure.  This  research  was  focused  on  the  study  of  the  fundamental 
processes  and  causes  of  frost  heaving  and  thaw-weakening  in  soils.  A  model  describing  frost  heaving  in  soils  has  been  developed  that  well 
captures  the  heaving  process  as  function  of  the  initial  and  boundary  conditions.  The  model  has  been  extended  to  include  the  entire 
freeze-thaw  cycle,  and  to  predict  the  consequences  of  phase  change  on  thaw-weakening  of  soils.  The  effects  in  soils  with  non-segregation 
freezing  have  been  included  in  the  model.  Frost  heaving  is  described  by  introducing  an  ice  growth  tensor,  which  captures  the  anisotropic 
growth  of  ice  lenses.  The  evolution  of  the  soil  strength  during  both  freezing  and  thawing  is  described  through  a  hardening/softening 
plasticity  model  including  the  effects  associated  with  freeze-thaw  cycles.  The  model  has  been  calibrated,  it  was  implemented  in  the  finite 
element  method,  and  it  was  used  to  simulate  boundary  value  problems. 

Enter  List  of  papers  submitted  or  published  that  acknowledge  ARO  support  from  the  start  of 
the  project  to  the  date  of  this  printing.  List  the  papers,  including  journal  references,  in  the 
following  categories: 

(a)  Papers  published  in  peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


01/27/2013  9.00  A.  DRESCHER,  R.L.  MICHALOWSKI.  Three-dimensional  stability  of  slopes  and  excavations, 

Geotechnique,  (01  2009):  0.  doi:  10.1680/geot.8.P.136 

01/27/2013  11.00  RadoslawL.  Michalowski.  Expanding  collapse  in  partially  submerged  granular  soil  slopes, 

Canadian  Geotechnical  Journal,  (12  2009):  0.  doi:  10.1 139/T09-064 

01/27/2013  10.00  Radoslaw  L.  Michalowski.  Limit  Analysis  and  Stability  Charts  for  3D  Slope  Failures, 

Journal  of  Geotechnical  and  Geoenvironmental  Engineering,  (04  2010):  583.  doi: 

08/11/201 1  2.00  Radoslaw  L.  Michalowski,  Tabetha  Martel.  Stability  Charts  for  3D  Failures  of  Steep  Slopes  Subjected  to 

Seismic  Excitation, 

Journal  of  Geotechnical  and  Geoenvironmental  Engineering,  (02  2011):  0.  doi:  10.1061/ 

(ASCE)GT.  1 943-5606.000041 2 

TOTAL:  4 


Number  of  Papers  published  in  peer-reviewed  journals: 


(b)  Papers  published  in  non-peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


TOTAL: 


Number  of  Papers  published  in  non  peer-reviewed  journals: 


(c)  Presentations 


Number  of  Presentations:  0.00 

Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


Received  Paper 


TOTAL: 


Number  of  Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 

Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


Received  Paper 


08/09/2011  4.00  Radoslaw  L.  Michalowski,  Srinivasa  S.  Nadukuru.  Delayed  increase  in  cone  penetration  resistance  of 

sand  after  dynamic  compaction, 

Second  International  Symposium  on  Computational  Geomechanics  (ComGeo  II),  Cavtat-Dubrovnik, 
Croatia..  2011/04/27  00:00:00,  .  :  , 

08/09/201 1  5.00  Srinivasa  S.  Nadukuru,  Tabetha  Martel,  Radoslaw  L.  Michalowski.  3D  analysis  of  steep  slopes  subjected 

to  seismic  excitation, 

Geo-Frontiers  2011,  Dallas,  TX..  201 1/03/13  00:00:00,  .  :  , 

08/09/201 1  6.00  Radoslaw  L.  Michalowski,  Srinivasa  S.  Nadukuru.  Stress  corrosion  cracking  and  delayed  increase  in 

penetration  resistance  after  dynamic  compaction  of  sand, 

Geo-Frontiers  201 1 ,  Dallas,  TX..  2011/03/14  00:00:00,  .  :  , 

08/09/201 1  7.00  Radoslaw  L.  Michalowski,  Srinivasa  S.  Nadukuru.  Post-liquefaction  state  of  sand,  stress  corrosion 

cracking,  and  relaxation  of  deviatoric  stress  in  previously  liquefied  sand  bed, 

5th  International  Conference  on  Earthquake  Geotechnical  Engineering  (5ICEGE).  Santiago,  Chile.. 

2011/01/10  00:00:00,.:, 

TOTAL:  4 


Number  of  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


(d)  Manuscripts 


Received 


01/27/2013  12.00  Yao  Zhang,  Radoslaw  L.  Michalowski.  Constitutive  model  and  simulation  of  non-segregation  freezing  and 

thawing  in  soils, 

16  Int.  Conf.  Soil  Mech.  Geotech.  Eng.  (12  2012) 

02/15/2013  13.00  Yao  Zhang,  Radoslaw  L.  Michalowski.  Thermal-mechanical  constitutive  modeling  for  freezing  and  thawing 

soils  , 

10th  International  Symposium  on  Cold  Regions  Development  (02  2013) 

08/09/201 1  3.00  Radoslaw  L.  Michalowski,  Yao  Zhang.  Frost-induced  heaving  of  soil  around  a  culvert, 

(08  2011) 

08/11/201 1  8.00  Radoslaw  L.  Michalowski,  Srinivasa  S.  Nadukuru.  Static  fatigue,  time  effects,  and  delayed  increase  in 

penetration  resistance  after  dynamic  compaction  of  sands, 

Journal  of  Geotechnical  and  Geoenvironmental  Engineering  (08  2011) 

TOTAL:  4 


Number  of  Manuscripts: 

Books 


Received  Paper 


TOTAL: 


Patents  Submitted 


Patents  Awarded 


Awards 

MTS  Visiting  Professor  of  Geomechanics,  University  of  Minnesota,  2012 
APEA  Tribute  Award  (American  Polish  Engineering  Association),  2010 

The  J.S.  Braun/Braun  Intertec  Professorship  in  Science  and  Technology,  University  of  Minnesota,  Twin  Cities,  2008 


Graduate  Students 


NAME 

PERCENT  SUPPORTED  Discipline 

Yao  Zhang 

0.75 

Srinivasa  Nadukuru 

0.00 

FTE  Equivalent: 

0.75 

Total  Number: 

2 

Names  of  Post  Doctorates 


NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Names  of  Faculty  Supported 


NAME 

PERCENT  SUPPORTED  National  Academy  Member 

Radoslaw  L.  Michalowski 

0.05 

FTE  Equivalent: 

0.05 

Total  Number: 

1 

Names  of  Under  Graduate  students  supported 


NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Student  Metrics 

This  section  only  applies  to  graduating  undergraduates  supported  by  this  agreement  in  this  reporting  period 


0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 


The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period: 
The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period  with  a  degree  in 

science,  mathematics,  engineering,  or  technology  fields:- 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  continue 
to  pursue  a  graduate  or  Ph.D.  degree  in  science,  mathematics,  engineering,  or  technology  fields: 

Number  of  graduating  undergraduates  who  achieved  a  3.5  GPA  to  4.0  (4.0  max  scale): . 
Number  of  graduating  undergraduates  funded  by  a  DoD  funded  Center  of  Excellence  grant  for 

Education,  Research  and  Engineering: . 
The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  intend  to 

work  for  the  Department  of  Defense 
The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  receive 
scholarships  or  fellowships  for  further  studies  in  science,  mathematics,  engineering  or  technology  fields: 


Names  of  Personnel  receiving  masters  degrees 

NAME 

Tabetha  Martel 

Total  Number:  1 

Names  of  personnel  receiving  PHDs 

NAME 

Total  Number: 


Names  of  other  research  staff 


NAME 


PERCENT  SUPPORTED 


FTE  Equivalent: 
Total  Number: 


Sub  Contractors  (DD882) 


Inventions  (DD882) 


See  Attachments 


Scientific  Progress 


Technology  Transfer 


Multi-Scale  Process  of  soil  Freezing,  Thawing, 
and  Thaw-Settlement 


Final  Report  to  U.S.  Army  Research  Office 
(Funding  W911NF-08- 1-0376,  Proposal  Number  53315-EV) 


Table  of  Contents 

1  Statement  of  the  problem  studied .  2 

2  Summary .  2 

3  Introduction  and  background . 2 

4  Development  of  the  constitutive  model  for  frost  heaving . 4 

4. 1  Porosity  rate  function  for  frost  heave  description .  4 

4.2  Simulation  of  the  frost  heave  around  a  culvert  using  PRF . 13 

4.3  Results  and  Discussions . 16 

5  Elasto-plastic  constitutive  model  for  freezing  and  thawing  soils . 19 

5.1  The  model . 19 

5.2  Thermal-mechanical  load  process . 23 

5.3  Implementation  and  application  of  the  model . 24 

6  Summary  of  the  most  important  results . 27 

References . 27 


Key  words: 

Freeze-thaw  cycle,  frost  heave,  ice  growth,  ice  lensing,  phase  change,  soil  thawing,  thaw¬ 
weakening,  mass  transfer,  heat  transfer,  frost  susceptibility,  constitutive  model 


List  of  Figures 


Figure  1.  Distribution  of  ice  fraction  (Test  DO . 7 

Figure  2.  Distribution  of  unfrozen  water  fraction  (Test  DO . 7 

Figure  3.  Retardation  factor  exp(-0i/0w)  (Test  DO . 8 

Figure  4.  Calibration  curves  for  clay  used  in  Fukuda  et  al.  (1997)  tests;  each  curve 

represents  the  simulation  for  a  different  thermal  gradient . 8 


Figure  5.  Comparison  of  the  total  frost  heave  simulated  by  the  model  developed,  one 

other  model,  and  the  experimental  tests  (for  different  temperature  gradients) .  9 


Figure  6.  Calibration  curves  for  Fukuda  et  al.  (1997)  tests  with  different  overburden 

pressure  (revised  model) . 10 

Figure  7.  Comparison  of  the  total  frost  heave  simulated  by  the  model  developed,  one 

other  model,  and  the  experimental  tests  (for  different  overburden  pressures)  10 

Figure  8.  JGST  freezing  process . 11 

Figure  9.  Test  data  and  calibration  curve  for  unfrozen  water  content  in  Alaskan  silt . 12 

Figure  10.  Calibration  of  the  model  with  JGST-freezing  test  for  Alaskan  silt . 13 

Figure  11.  Geometry  and  mesh  of  the  model  culvert . 13 

Figure  12.  The  annual  air  temperature  variation  in  Aniak,  AK,  1961  -  1962 .  14 

Figure  13.  (a)  Below-ground  temperature  in  Aniak,  AK,  1952  -  1953  (after  Aitken  and 

Ful wider,  1962),  (b)  Initial  temperature  profile  for  ID  simulation . 15 

Figure  14.  Temperature  regimes  for  ID  simulations . 15 

Figure  15.  Frost  heave  simulation  results  (1-D) . 15 

Figure  16.  Initial  temperature  distribution  at  t  =  0 . 16 

Figure  17.  (a)  Temperature  distribution  after  37  days,  and  (b)  frost  heave  contour . 17 

Figure  18.  (a)  Temperature  distribution  after  100  days,  and  (b)  frost  heave  contour . 17 


Figure  19.  (a)  Temperature  distribution  after  46  days,  and  (b)  frost  heave  contour  (culvert 
with  outside  thermal  insulation  around  the  circumference,  X  =  0.02  W/(m°C) ). 
. 18 

Figure  20.  (a)  Temperature  distribution  after  100  days,  and  (b)  frost  heave  contour 


(culvert  with  outside  thermal  insulation  around  the  circumference,  X  =  0.02 
W/(m°C) ) . 18 

Figure  21.  (a)  Freeze-thaw  thermal  loading  path  in  v-ln (p)  plane;  (b)  “pseudo 

preconsolidation  pressure”  changing  along  BD  due  to  increasing  pore  ice  ratio 
. 21 

Figure  22.  Freeze-thaw  cycle  and  load  path  with  corresponding  yield  curves . 24 

Figure  23.  Calibration  of  preconsolidation  stress  for  frozen  soil  vs  temperatures . 25 

Figure  24.  Boundary  condition  on  top  of  the  soil  column . 26 

Figure  25.  Temperature  distribution  at  t  =  0  and  t  =  300  h . 26 

Figure  26.  Compression  for  element  #18  (a)  and  element  #4  (b) . 27 


II  P  a  g  e 


1  Statement  of  the  problem  studied 

Phenomena  associated  with  freeze-thaw  cycles  in  soils  are  not  fully  understood,  which 
hinders  the  development  of  efficient  models  for  model-based  simulations.  Both  the 
formation  of  ice  lenses  during  freezing  and  soil  weakening  during  thawing  were  studied. 
In  particular,  a  porosity  growth  function  was  considered  with  the  goal  of  developing  a 
frost  heave  model,  and  a  plasticity-based  solid  model  with  hardening  and  softening 
governed  by  ice  content  was  studied  to  model  increase  in  strength  during  freezing  and 
weakening  during  thawing. 

2  Summary 

Thermal  effects  in  soils  are  relevant  to  both  Army  and  civilian  operations,  particularly  in 
the  regions  of  seasonal  freezing  and  permafrost.  The  significance  of  thermal  processes  in 
soils  may  be  amplified  if  the  current  trends  in  climate  changes  continue.  Of  problems  that 
require  attention  are  the  frost  heave  and  thaw  weakening  of  soils  during  freeze-thaw 
cycles.  Both  have  been  studied  with  the  goal  of  developing  a  mathematical  description 
that  could  be  used  successfully  in  model-based  simulations. 

The  progress  in  the  description  of  frost  heave  has  been  hindered  by  the  lack  of 
understanding  of  the  factors  behind  development  of  cryogenic  suction  in  freezing  soils. 
From  the  numerical  standpoint,  on  the  other  hand,  the  difficulties  lay  in  the  fact  that  some 
of  the  crucial  parameters  in  the  freezing  zone  evolve  by  orders  of  magnitude  (e.g., 
hydraulic  conductivity),  and  evaluation  (calibration)  of  those  parameters  are  quite 
difficult.  In  modeling,  this  problem  was  overcome  by  introducing  a  concept  of  porosity 
rate  function,  which  allows  eliminating  hydraulic  conductivity  from  the  model  by 
introducing  a  porosity  growth  tensor,  which  was  successfully  calibrated. 

The  second,  equally  important  issue  is  that  of  the  evolution  of  strength  of  the  soil  during 
freezing  and  thawing.  A  model  was  developed  based  on  hardening  plasticity,  with  ice 
content  and  specific  volume  being  the  principal  parameters  responsible  for  hardening  and 
softening  of  the  soil.  The  models  were  used  with  success  in  simulations  of  freezing  and 
thawing  of  frost-susceptible  soils. 

3  Introduction  and  background 

Cold  regions,  including  permafrost  and  seasonal  freezing  areas  constitute  a  large  portion 
of  Earth’s  surface.  In  the  northern  hemisphere,  more  than  20%  of  land  surface  is  occupied 
by  permafrost,  whereas  more  than  50%  experiences  seasonal  freezing  and  thawing.  Frost 
action  in  soils  involving  frost  heaving  and  thaw  settlement  is  commonly  seen  in  seasonal 
freezing  areas  rendering  built  infrastructure  vulnerable  to  damage.  Broken 
communication  channels,  damaged  lifelines,  broken  pipelines,  malfunctioning  utilities, 
cracked  pavements,  jacked  up  bridge  foundations,  tilted  structures,  are  all  examples  of 
damage  suffered  from  frost  action.  Long-term  paleographic  records  indicate  an  on-going 
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warming  of  the  climate,  which  has  resulted  in  thawing  of  portions  of  the  permafrost  area. 
The  increase  in  permafrost  temperature  leads  to  thickening  of  the  active  layer  (upper 
crust  layer  where  active  freezing-thawing  cycles  occur),  leading  to  extensive  settlement 
of  the  ground  surface  causing  damage  to  infrastructure.  Understanding  how  the  soil 
behaves  upon  freezing  and  thawing  has  the  potential  of  changing  the  operation  practices 
and  design  philosophies,  and  developing  methods  to  alleviate  the  damages. 

Various  soils  experiencing  freezing  and  thawing  may  exhibit  dramatically  different 
behaviors.  A  good  portion  of  silty  soils  and  clays  are  characterized  by  susceptibility  to 
frost  heaving.  This  process  is  caused  by  transfer  of  moisture  from  unfrozen  layers  of  soil 
into  the  freezing  zone,  driven  by  gradient  in  cryogenic  suction,  and  freezing  to  form 
segregated  ice  referred  to  as  ice  lenses.  The  zone  where  this  process  occurs  is  often  called 
the  frozen  fringe.  The  frost  heave  of  a  frost-susceptible  soil  could  be  more  than  20%  of 
the  thickness  of  the  frozen  layer,  and  it  depends  on  the  rate  of  freezing  front  propagation 
through  this  layer.  Upon  thawing,  the  thaw  settlement  could  be  larger  or  smaller  than  the 
displacement  caused  by  frost  heave,  depending  on  whether  the  soil  was  normally 
consolidated  ( NC)  or  over-consolidated  (OC)  prior  to  freezing. 

Thawing  of  a  frozen  NC  soil  in  its  first  freeze-thaw  cycle  could  cause  a  settlement  larger 
than  the  heave  induced  during  the  freezing  phase.  However,  some  heave  could  remain  in 
an  OC  soil  subjected  to  its  first  freeze-thaw  cycle.  Artificial  ground  freezing  applied  in 
soft  soil  construction  (tunneling,  excavations),  and  pipelines  transporting  chilled  media 
through  unfrozen  soil,  are  examples  of  circumstances  where  soil  may  be  subjected  to  its 
first  freeze-thaw  cycle. 

For  a  non-frost-susceptible  soil,  such  as  sand  or  gravel,  no  ice  segregation  will  take  place 
during  freezing.  Pore  ice  will  be  formed  due  to  the  phase  change  of  the  water  in  the  pores. 
While  the  transfer  of  moisture  can  still  occur  due  to  the  thermal  gradient,  ice  lenses  do 
not  form,  and  the  macroscopic  deformation  of  the  soil  mixture  (mineral,  ice,  water,  air) 
upon  freezing  is  owed  to  expansion  of  freezing  water.  This  expansion  may  be  fully 
accommodated  by  pores  in  unsaturated  soils,  but  it  will  cause  macroscopic  deformation 
in  water-saturated  soils 

In  addition  to  the  deformation  induced  by  the  freezing  process  in  soils,  the  strength  gain 
(hardening)  during  freezing  and  soil  weakening  (softening)  during  thawing,  are  important 
concerns  in  cold  regions  operations.  Experimental  data  indicates  that  an  increase  in 
strength  occurs  in  both  the  frost-susceptible  and  non-frost-susceptible  soils  as  they  freeze. 
A  reduction  in  strength  is  then  expected  upon  thawing. 

Constitutive  models  describing  changes  in  strength  caused  by  freeze-thaw  cycles  have 
not  been  developed,  though  some  efforts  have  been  published  recently  (e.g.,  Lai  et  al., 
2009,  Wei  et  al.,  2011,  Shastri  and  Sanchez,  2012).  These  models  assume  elasto-plastic 
behavior  with  account  for  viscous  properties  of  ice  (creep);  they  are  based  on  continuum 
approach,  and  are  customized  for  one  type  of  the  frozen  soil.  Changes  of  the  soil 
components  upon  freezing  and  thawing,  and  the  corresponding  changes  in  strength,  are 
not  addressed  by  these  models. 
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The  model  developed  in  this  research  contains  two  parts:  a  revised  porosity  rate  model 
for  frost  heaving  and  an  elasto-plastic  constitutive  model  with  evolution  of  soil  strength 
caused  by  freezing  and  thawing.  The  constitutive  model  developed  is  capable  of 
addressing  freezing  and  thawing  effects,  it  captures  the  deformation  behavior,  and 
describes  the  evolution  of  strength  in  the  soil  subjected  to  loading  and  thermal  changes, 
including  phase  change.  The  model  is  based  on  the  critical  state  concept  and  it  is 
formulated  by  introducing  the  influence  of  ice  ratio  into  the  hardening/softening 
plasticity  framework.  The  model  was  implemented  in  the  finite  element  method,  and 
boundary  value  problems  were  successfully  simulated. 

4  Development  of  constitutive  model  for  frost  heaving 

4.1  Porosity  rate  function  for  frost  heave  description 

Frost  heave  is  caused  by  the  growth  of  ice  lenses,  and  the  porosity  rate  function  describes 
the  average  growth  of  ice  in  a  soil  volume.  This  model  was  found  to  be  very  effective  in 
simulations  of  frost  heave  processes.  A  great  advantage  of  this  model  is  a  relatively  small 
number  of  parameters  needed  to  describe  the  process  of  frost  heave.  The  model  assumes 
the  soil  to  be  saturated,  and  the  soil  skeleton  to  behave  as  an  elastic  solid.  An  early 
proposal  of  a  porosity  rate  function,  central  to  this  model,  was  suggested  in  Michalowski 
(1993),  and  it  was  modified  later  by  Michalowski  and  Zhu  (2006).  The  porosity  growth 
in  this  model  is  described  by  the  porosity  rate  tensor 

=  '»  4  (1) 

where  st  is  a  unit  strain  tensor  that  maps  the  ice  volume  increase  into  anisotropic  growth 

consistent  with  the  growth  of  ice  lenses,  and  h  is  a  scalar  function  representing  porosity 
growth 
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Notation:  h  is  the  porosity  rate  ( dn/dt ),  T  is  the  temperature  of  the  soil  (°C),  and  au.  is 
the  first  invariant  of  the  stress  tensor  in  the  frozen  soil  (zero  if  tension).  Parameters  hm 
and  Tm  are  the  maximum  porosity  rate  and  the  temperature  (°Q  at  which  that  maximum 
porosity  rate  occurs,  respectively;  T0  is  the  freezing  temperature  of  water;  hm  is 
determined  at  well-defined  temperature  gradient  gT;  £  is  a  material  parameter  defining 


stress  dependency.  The  quotient 


ST/ 


c'l 


/  gT  indicates  linear  relationship  of  the  porosity 


rate  to  the  temperature  gradient  in  the  heat  flow  direction.  The  term  exp(-|cr/d |/ <^)is  a 
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retardation  function:  frost  heave  can  be  significantly  reduced  by  the  compressive  stress  in 
the  direction  of  heat  flow,  with  some  contribution  of  the  normal  stresses  in  the  other  two 
perpendicular  directions  (parameter  C,  was  calibrated  for  specific  clay).  A  porosity 
threshold  of  0.75  was  introduced,  past  which  further  heave  ceases,  to  indicate  the  fact  that 
the  frost  heave  stops,  or  reduces  to  an  insignificant  rate,  after  the  freezing  front  reaches  a 
steady-state  location  (Fukuda  et  al.,  1997). 

A  modification  of  this  model  was  proposed  to  make  it  more  realistic  in  describing  the 
influence  of  stress  state  and  the  soil  components  (soil  skeleton,  unfrozen  water  and  ice)  to 
the  frost  heave  process.  The  modified  porosity  rate  function  has  the  following  form 
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where  a \  is  the  stress  component  in  the  heat  flow  direction,  and  0.  and  0W  are  the 
volumetric  fractions  of  ice  and  unfrozen  water,  respectively. 

The  first  invariant  of  the  total  stress  tensor  in  soil  akk  was  replaced  by  the  stress  in  the 
heat  flow  direction  a) .  The  reason  for  that  modification  is  that  the  ice  lenses  tend  to  grow 

in  the  heat  flow  direction  and  it  is  the  stress  component  in  that  direction  that  inhibits  the 
heave. 

In  frozen  silt  and  clay,  only  a  portion  of  liquid  water  turns  into  solid  phase  at  the  freezing 
point,  and  the  rest  remains  in  a  liquid  state  at  the  temperature  well  below  freezing.  A 
3-parameter  function  (Michalowski  1993)  is  chosen  here  to  describe  the  unfrozen  water 
content  in  the  frozen  soil 


w  =  w*  +  (w  —  w*)  ea(T  To)  (4) 

The  calibration  for  clay  in  Fukuda  et  al.  (1997)  tests  results  in  the  following  set  of 
parameters:  w*  =  0.058,  w  =  0.285,  T0  =  0 °C,  and  a  =  1.6.  In  the  following,  the 

modified  porosity  rate  function  is  calibrated  for  clay  based  on  an  extensive  set  of  tests  by 
Fukuda  et  al.  (1997). 

The  clay  tested  by  Fukuda  et  al.  was  saturated,  with  the  initial  porosity  of  0.427,  and 
specific  gravity  of  2.62.  The  density  of  the  soil  was  as  1,927  kg/m3,  the  Young’s  modulus 
and  Poisson’s  ratio  were  taken  as  E  =  1 1.2  MPa  and  //  =  0.3  for  both  frozen  and  unfrozen 
soil.  The  typical  values  of  thermal  properties  for  soil  particles,  water  and  ice  are  listed  in 
Table  1.  These  values  are  used  in  both  numerical  calibration  tests  and  in  the  numerical 
simulations  of  boundary  value  problems. 
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Table  1  Thermal  properties  of  clay  soil  constituents  (after  Williams  and  Smith,  1989) 


Density 

P(kg/m3) 

Mass  heat 

Volumetric  heat 

Latent 

Thermal 

capacity 

capacity 

heat 

Conductivity 

c(J/(kg-°0) 

C  (J/(m3-°C)) 

L  (J/kg) 

X  (W/m3-°C) 

Soil  particles 
(clay  mineral) 

2620 

900 

2.36T06 

2.92 

Water 

1000 

4180 

4.18-106 

3.33-105 

0.56 

Ice 

917 

2100 

1.93-106 

2.24 

The  porosity  rate  function  in  equation  (3)  was  implemented  in  the  FEM  system 
ABAQUS,  with  a  porosity  growth  tensor  and  thermal  conservation  equations  which  can 
be  found  in  Michalowski  and  Zhu  (2006).  Then  the  implemented  model  was  used  to 
simulate  a  set  of  soil  freezing  tests  performed  by  Fukuda  et  al.  (1997).  The  specimens 
were  cylindrical,  100  mm  in  diameter  with  initial  height  of  70  mm.  The  initial  and 
boundary  conditions  for  the  tests  performed  under  constant  load  are  listed  in  Table  2.  All 
these  tests  had  ramped  temperatures.  The  frost  heave  parameters  for  the  particular  clay 
used  in  Fukuda  et  al.  (1997)  tests  were  found  to  be:  Tm  =  -0.74°C,  nm  =  1.88TO'5  s'1,  gT 

=  100°  C/m,  and  £  =  0.42  MPa. 


Table  2  Boundary/initial  conditions  for  tests  by  Fukuda  et  al.  (1997) 


Ramped  Freezing  Tests 
(dT/dl) 

Warm  plate  (top) 
(°C) 

Cold  plate  (bottom) 
(°C) 

Overburden  pressure 
(kPa) 

A 

0.095 

7-0.042t* 

-0.042t 

25 

B 

0.070 

5-0.042t 

-0.042t 

25 

C 

0.061 

4-0.042t 

-0.042t 

25 

Di 

25 

d2 

150 

d3 

0.045 

3-0.042t 

-0.042t 

300 

d4 

400 

d5 

600 

E 

0.035 

2-0.042t 

-0.042t 

25 

- - - 

t-time  (in  hours),  total  testing  time  47hrs. 


_ 

To  illustrate  the  role  term  e  B"'  plays  in  the  porosity  rate  function,  the  volumetric  content 
of  ice,  water  and  the  value  of  the  retardation  coefficient  in  the  simulated  test  Di  (see 
Table  2)  are  plotted  in  Figure  1  through  Figure  3.  As  the  ice  fraction  increases,  the  term 

_  A 

e  e"  decreases.  This  is  the  effect  that  allowed  elimination  of  shut-off  porosity  threshold 
from  the  earlier  model. 
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Figure  1.  Distribution  of  ice  fraction  (Test  Di) 


Figure  2.  Distribution  of  unfrozen  water  fraction  (Test  Di) 
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Figure  3.  Retardation  factor  exp(-0i/0w)  (Test  DO 


The  calibration  (curve  fitting)  process  was  performed  first  using  tests  with  different 
temperature  gradients  (A,  B,  C,  Di,  E  in  Table  2),  and  then  the  influence  of  stress  was 
calibrated  using  the  tests  with  different  overburden  load  (Di  to  D5  in  Table  2).  The  first 
family  of  the  calibration  curves  (dependence  on  the  thermal  gradient)  is  shown  in  Figure 
4. 


Figure  4.  Calibration  curves  for  clay  used  in  Fukuda  et  al.  (1997)  tests;  each  curve 
represents  the  simulation  for  a  different  thermal  gradient 
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Comparison  of  the  heave  magnitudes  from  the  best  calibration  (curve  fitting)  effort  is 
illustrated  in  Figure  5. 


Figure  5.  Comparison  of  the  total  frost  heave  simulated  by  the  model  developed,  one 
other  model,  and  the  experimental  tests  (for  different  temperature  gradients) 


The  respective  graphs  from  the  calibration  effort  with  the  tests  of  varied  overburden 
pressure  are  shown  in  Figure  6  and  Figure  7. 
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15 


Figure  6.  Calibration  curves  for  Fukuda  et  al.  (1997)  tests  with  different  overburden 

pressure  (revised  model) 


Figure  7.  Comparison  of  the  total  frost  heave  simulated  by  the  model  developed,  one 
other  model,  and  the  experimental  tests  (for  different  overburden  pressures) 


The  Root-Mean-Square  Error  (RMSE)  is  0.0098  to  0.0206  for  calibration  with  different 
temperature  gradients,  and  0.0052  to  0.0069  for  calibration  with  respect  to  different 
overburden  load. 
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Additional  calibration  was  performed  for  Alaskan  silt,  based  on  tests  performed 
according  to  Japanese  Geotechnical  Standard  Test  (JGST)  by  Kim  (2011)  at  the 
Univsersity  of  Alaska,  Fairbanks.  JGST  freezing  test  is  the  standard  test  method  for 
obtaining  frost  sussceptibility  in  Japan  (Japan  Geotechnical  Society  2003),  which 
requires  the  temperature  of  the  warm  boundary  to  remain  constant  slightly  above  0°C, 
whereas  the  temperature  of  the  cold  boundary  is  gradually  reduced  at  a  constant  rate. 
Figure  8  illustrates  the  temperature  change  for  a  JGST  freezing  test. 


Figure  8.  JGST  freezing  process 

The  set  of  tests  selected  to  calibrate  the  model  are  JGST9  to  JGST12  for  Alaska  silts 
(Kim,  2011).  The  test  conditions  are  summarized  in  Table  3.  Since  the  set  of  JGST  frost 
heave  tests  selected  contains  effects  of  both  different  overburden  pressure  and 
temperature  gradient,  parameters  in  porosity  rate  function  can  be  determined. 

Table  3  Conditions  for  JGST  tests 


Test  No. 

Linear  reduction  of  top 
pedestal  temperature  (°C) 

Average  bottom 
pedestal  temperature 

ra 

Operation 

time 

(h) 

Overburden 

pressure 

(kPa) 

JGST-  9 

-0.26  ->  -3.93 

0.31 

30 

JGST-10 

-0.19  ->  -3.80 

0.16 

40 

40 

JGST-11 

-0.26  -3.86 

0.32 

60 

JGST-12 

-0.24  ->  -3.84 

0.29 

80 

The  calibration  for  unfrozen  water  content  in  Alaskan  silt  results  in  the  following  set  of 
parameters:  w*  =  0.813,  w  =  0.325,  T0  =  0°C,  and  a  =  6.0. 
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Figure  9.  Test  data  and  calibration  curve  for  unfrozen  water  content  in  Alaskan  silt 

The  Alaskan  silt  was  saturated,  with  the  initial  porosity  of  0.47,  and  specific  gravity  of 
2.62.  The  density  of  the  soil  was  1,923  kg/m3,  the  Young’s  modulus  and  Poisson’s  ratio 
were  taken  as  E  =  1 1.2  MPa  and  n  =  0.3  for  both  frozen  and  unfrozen  soil.  The  properties 
for  Alaskan  silt  particles  were  measured  by  Kim  (2011)  and  are  listed  in  Table  4. 
Properties  for  ice  and  water  are  given  in  Table  1. 

Table  4  Thermal  properties  of  Alaska  silt  soil  particles  (Kim,  2011) 


Density 

p(kg/m3) 

Mass  heat 

Volumetric  heat 

Thermal 

capacity 

capacity 

Conductivity 

c(J/(kg-°C)) 

C  (J/(m3-°C)) 

A  (W/m3-°C) 

Soil  particles 
(clay  mineral) 

2746 

800 

2.20-106 

2.92 

The  calibration  with  Kim  (2011)  is  illustrated  in  Figure  10.  The  parameters  obtained  from 
calibration  are:  Tm  =  -0.20°C,  hm  =  1.80T0"5  s'1,  gT  =  100°C/m,  and  £  =  0.058  MPa. 

These  parameters  will  be  used  in  a  simulation  to  illustrate  the  application  of  the  model, 
and  to  reveal  an  interesting  phenomenon  of  a  frost-induced  heave  around  a  culvert. 
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Figure  10.  Calibration  of  the  model  with  JGST-freezing  test  for  Alaskan  silt 


4.2  Simulation  of  the  frost  heave  around  a  culvert  using  PRF 

To  illustrate  applicability  of  the  model,  simulations  of  freezing  of  an  unpaved  road  above 
a  culvert  in  frost  susceptible  soil  were  carried  out.  The  model  was  implemented  in  the 
FEM  system  ABAQUS.  The  geometry  of  the  boundary  value  problem  and  the  finite 
element  mesh  of  the  culvert  and  the  surrounding  soil  are  shown  in  Figure  11.  The  top  of 
the  culvert  in  this  simulation  is  located  1.0  m  below  the  ground  surface.  The  external 
diameter  of  it  is  1.0  m  and  its  wall  thickness  is  10.0  mm.  The  material  of  the  culvert  is 
steel,  and  a  perfect  bonding  with  soil  (ad-frozen  interface)  is  modeled.  The  simulations 
were  performed  for  the  boundary  conditions  consistent  with  Alaska  climate  at  Aniak,  AK. 


Figure  1 1 .  Geometry  and  mesh  of  the  model  culvert 
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More  details  on  this  simulation  can  be  found  in  Zhang  and  Michalowski  (2012).  The  data 
of  daily  air  temperature  change  in  Aniak,  AK,  from  year  1961  to  1962  was  obtained  from 
National  Climate  Data  Center  and  it  is  shown  in  Figure  12.  The  trend  of  the  seasonal 
temperature  change  is  very  clear  and  regular.  The  highest  temperature  is  around  20°C  and 
the  lowest  temperature  is  around  -40°C.  Freezing  temperatures  last  for  about  2/3  of  a  year. 
When  freezing  starts,  the  temperature  drops  down  from  above  zero  to  about  -30°C  in 
about  3  months.  The  mean  daily  temperature  amplitude  (defined  as  the  difference 
between  the  peak  and  the  mean  daily  value)  is  about  10°C,  although  it  happens 
sometimes  to  reach  as  much  as  20°C. 
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Figure  12.  The  annual  air  temperature  variation  in  Aniak,  AK,  1961  -  1962 
(National  Climate  Data  Center) 

The  ground  temperature  in  Alaska  was  monitored  by  U.S.  Army  Cold  Regions  Research 
and  Engineering  Laboratory  (CRREL)  researchers  during  the  year  1947  to  1958  at 
several  locations.  The  ground  temperature  in  Aniak,  AK,  was  reported  and  the 
temperature  profile  below  the  ground  surface  in  shown  in  Figure  13(a)  from  the  end  of 
freezing  season  in  1952  until  the  middle  of  the  freezing  season  in  1953  (Aitken  and 
Fulwider,  1962).  This  figure  indicates  small  variations  of  temperature  at  about  4  m  (12  ft) 
below  the  ground  surface,  and  the  temperature  is  nearly  constant  at  about  3°C  (37°F)  at 
depth  of  about  7  m  (22  ft). 
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(a)  GROUND  TEMPERATURE  IN  DEGREES  CENTIGRADE 
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Figure  13.  (a)  Below-ground  temperature  in  Aniak,  AK,  1952  -  1953  (after  Aitken  and 
Fulwider,  1962),  (b)  Initial  temperature  profile  for  ID  simulation 

A  portion  of  the  temperature  function  from  early  October,  1952  (day  270)  to  early 
January,  1953  (day  370),  a  total  100  days  (shown  in  Figure  12)  were  used  in  the 
simulation.  The  bottom  boundary  of  the  model  in  Figure  1 1  is  assumed  to  have  a  constant 
temperature  of  3°C.  The  thermal  boundaries  between  air  and  soil,  air  and  culvert  are 
simulated  using  a  Fourier  boundary  condition 


XVT  =  ht(Tmr-T) 


(5) 


where  X  and  hc  are  the  soil  heat  conductivity  and  convective  heat  transfer  coefficient,  T  is 
soil  temperature  and  Tair  is  the  ambient  temperature. 


Figure  14.  Temperature  regimes  for  ID 
simulations 


Figure  15.  Frost  heave  simulation  results 
(1-D) 


One-dimensional  soil  freezing  tests  were  simulated  first,  using  the  true  temperature  data 
and  the  linearized  temperature  change  (curve  1  and  2  in  Figure  14).  These  temperature 
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functions  are  used  as  the  boundary  conditions  at  the  top  of  a  seven-meter  column  in  ID 
simulation,  whereas  the  bottom  temperature  is  kept  at  3°C.  The  initial  temperature 
distribution  for  the  ID  simulation  is  shown  in  Figure  13(b).  The  ID  frost  heave 
simulation  results  are  shown  in  Figure  15.  The  results  indicate  that  the  ID  frost  heave 
simulated  using  the  linearized  temperature  function  is  quantitatively  similar  to  that  from 
simulation  with  the  true  temperature. 
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Figure  16.  Initial  temperature  distribution  at  t  =  0. 

The  linearized  air  temperature  function  was  used  in 
the  culvert  simulation.  At  time  t  =  0,  the 
temperature  of  the  ground  surface  and  the 
temperature  inside  the  culvert  are  set  to  be  5°C. 
Then  the  ambient  temperature  drops  down  to  -30°C 
in  100  days  (shown  in  Figure  15  line  2).  The 
temperature  at  the  base  of  the  model  (bottom 
boundary)  is  kept  constant  at  3°C.  The  initial 
temperature  distribution  is  shown  in  Figure  16. 


The  values  of  the  parameters  in  the  porosity  rate  function  used  here  were  taken  from 
calibration  based  on  tests  done  by  Kim  (2011)  and  they  were  given  in  the  earlier  Section 
describing  the  porosity  rate  function.  However,  to  simulate  the  soil  that  is  less  frost 
susceptible,  the  value  of  nmwas  set  to  be  10%  of  the  value  from  the  calibration:  nm  = 

1.80T0'6  s4 


The  soil  is  assumed  to  be  saturated,  with  the  initial  porosity  of  0.47,  and  specific  gravity 
of  2.62.  The  density  of  the  soil  was  1,923  kg/m3,  the  Young’s  modulus  and  Poisson’s 
ratio  were  taken  as  £=11.2  MPa  and  //  =  0.3  for  both  frozen  and  unfrozen  soil.  The 
thermal  parameters  were  taken  from  Williams  and  Smith  (1989)  and  Selvadurai  et  al. 
(1999).  The  thermal  conductivities  for  soil  particles,  water,  and  ice  were  taken  as  2.92, 
0.56,  and  2.24  W/(m°C).  The  heat  capacities  for  soil  particles,  water,  and  ice  were  set  to 
800,  4180,  and  2100  J/(kg°C).  The  latent  heat  of  fusion  of  water  is  3.33-105  (J/kg).  The 
properties  of  the  culvert  steel  were  taken  as:  thermal  conductivity  36  W/(m°C),  mass  heat 
capacity  477.3  J/(kg°C),  and  density  7,800  kg/m3. 

In  addition  to  simulating  the  boundary  value  problem  defined  in  Figure  11,  an  influence 
of  using  culvert  insulation  on  the  frost  heave  resulting  from  the  same  thermal  boundary 
conditions  was  studied.  A  10  mm  thick  layer  of  insulation  with  a  heat  conductivity  of  X  = 
0.02  W/(m°C)  was  placed  around  the  outside  perimeter  of  the  culvert.  The  density  of  the 
insulation  is  50  kg/m3,  the  Young’s  modulus  and  Poisson’s  ratio  were  taken  as  E  =  10.0 
MPa  and  ju  =  0.3,  and  the  heat  capacity  was  2,000  J/(kg°C). 

4.3  Results  and  discussions 


The  results  of  frost  heave  simulation  around  a  culvert  are  shown  in  Figure  17  and  Figure 
18.  At  first,  a  “bump”  is  formed  at  the  surface  directly  above  the  culvert;  this  bump 
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reaches  its  maximum  after  37  days.  Thereafter,  the  bump  stops  growing  and  turns  into  a 
“dip”;  the  displacements  after  100  days  are  shown  in  Figure  18. 
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Figure  17.  (a)  Temperature  distribution  after  37  days,  and  (b)  frost  heave  contour. 
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Figure  18.  (a)  Temperature  distribution  after  100  days,  and  (b)  frost  heave  contour. 

When  the  air  temperature  drops  below  freezing,  two  freezing  fronts  propagate  into  the 
soil  directly  above  the  culvert:  one  down  from  the  surface,  and  the  second  one  up  from 
the  culvert.  If  the  soil  is  frost-susceptible,  the  heave  process  (ice  lens  growth)  will  occur 
in  the  soil.  Initially,  this  process  is  particularly  intense  above  the  culvert,  as  this  portion 
of  the  soil  has  the  freezing  front  propagating  both  from  the  surface  of  the  soil  and  from 
the  culvert  side.  The  resulting  outcome  is  a  bump  on  the  surface.  If  the  freezing 
temperature  continues,  the  bump  will  continue  to  grow  until  all  soil  above  the  culvert 
becomes  all  frozen.  After  that,  the  heave  of  the  area  above  the  culvert  slows  down  (or 
drops  down  to  insignificant  rate),  but  the  heave  in  regions  further  away  from  the  culvert 
continues,  which  eventually  results  in  a  “dip”  in  the  surface  above  the  culvert. 
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The  key  factor  that  determines  the  profile  of  the  surface  above  the  culvert  is  the  boundary 
condition.  If  the  freezing  temperature  continues  only  for  a  short  period,  a  “bump”  will 
form  on  the  ground  surface,  whereas  if  the  freezing  temperature  persists  for  a  long  time,  a 
“depression”  is  likely  to  occur.  While  the  former  is  confirmed  by  observations 
(Andersland  and  Ladanyi  2004)  the  latter  is  a  hypothetical  long-term  outcome. 
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Figure  19.  (a)  Temperature  distribution  after  46  days,  and  (b)  frost  heave  contour  (culvert 
with  outside  thermal  insulation  around  the  circumference,  A  =  0.02  W/(m°C) ). 
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Figure  20.  (a)  Temperature  distribution  after  100  days,  and  (b)  frost  heave  contour 
(culvert  with  outside  thermal  insulation  around  the  circumference,  A  =  0.02  W/(m°C) ). 

The  results  of  frost  heave  simulation  around  culvert  with  thermal  insulation  are  shown  in 
Figure  19  and  Figure  20.  In  Figure  19,  the  maximum  “bump”  is  formed  at  about  46  days, 
and  the  frost  heave  is  smaller  compared  to  that  in  Figure  17.  However,  the  “depression” 
in  Figure  20  is  as  large  as  that  in  Figure  18. 
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5  Elasto-plastic  constitutive  model  for  freezing  and 
thawing  soils 

5.1  The  model 

Little  effort  was  devoted  in  the  past  to  constitutive  modeling  of  soils  subjected  to  freezing 
and  thawing,  and  no  reliable  models  have  been  developed,  so  far.  After  1980’s  some 
effort  was  made  toward  describing  reversible  and  irreversible  deformation,  but  these 
efforts  did  not  give  rise  to  computational  tools  that  could  be  used  in  model-based 
simulations.  Most  of  the  models  were  based  on  continuum  approach,  and  were  suited  for 
one  type  of  frozen  soil.  The  changes  in  strength  due  to  varying  soil  composition  upon 
freezing  and  thawing  were  not  addressed  by  these  models.  The  most  recent  efforts  draw 
from  the  area  of  unsaturated  soil  mechanics,  and  they  exploit  an  analogy  of  cryogenic 
suction  in  freezing  soil  to  matric  suction  in  unsaturated  soils  (e.g.,  Shastri  and  Sanchez, 
2012;  Nishimura,  et  al.,  2009).  The  development  of  a  model  that  accounts  for 
deformation  and  strength  evolution  caused  by  freezing  and  thawing  was  part  of  this 
project,  and  it  was  motivated  by  the  absence  of  models  that  could  be  effectively  used  in 
model-based  simulations 

The  model  developed  introduces  ice  ratio  ei  as  a  key  parameter  in  the  description  of  the 
behavior  of  the  frozen  soil.  The  ice  ratio  is  defined  as 


where  V,  is  the  volume  of  ice  and  Vs  is  the  volume  of  the  solid  constituent  (mineral 
skeleton)  in  a  given  volume  of  soil.  Ice  ratio  ej  is  uniquely  related  to  the  unfrozen  water 
content  (see  equation  (4)) 


ei  =  ( vv’o  - vv) '  ~  1  -09  (7) 

Pw 

where  wo  is  the  moisture  content  of  saturated  soil  before  freezing,  w  is  the  current 
unfrozen  moisture  content,  and  ps  and  p„.  are  the  density  of  the  mineral  of  the  skeleton 
and  water,  respectively.  Specific  volume  v  and  mean  stress  p  are  introduced 

V  1 

v  =  —  =  1  +  e’  P  =  ^akk  k  =  1,2,3  (8) 

where  V  is  the  total  volume  of  the  soil,  e  is  the  void  ratio,  and  a«-  is  the  first  invariant  of 
the  stress  state.  Compression  tests  on  frozen  soils  indicate  that  the  behavior  can  be 
represented  by  the  normal  compression  line  (NCL)  in  the  u,  In p  plane,  and  the  unloading¬ 
reloading  line  (URL)  (Qi  et  al.,  2010,  Lee  et  al.,  2002).  The  slopes  of  these  lines  vary, 
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and  they  depend  on  the  extent  of  freezing.  The  slopes  for  the  two  lines  are  defined  by 
coefficients  A  and  k  for  unfrozen  soil,  and  Af  and  Kf  for  frozen  soil,  respectively.  The 

specific  volume  upon  isotropic  compression  of  frozen  soil  is  given  by 

v  =  vf  -  Af  In  p  (9) 

and  the  elastic  behavior  in  unloading-reloading  regime  is  given  by 

dV  =  -k  ^  (10) 

P 

In  the  model  developed,  both  Af  and  Kf  are  functions  of  ice  ratio  The  yield  condition 

of  frozen  soil  is  described  here  with  an  ellipse  in  the  deviatoric-mean  stress  plane  ( q ,  p), 
and  its  analytical  expression  is  presented  later  in  equation  (19).  Such  a  surface  has  been 
in  use  for  quite  some  time  in  the  framework  of  the  critical  state  soil  mechanics,  and 
because  of  its  effectiveness  and  relative  simplicity,  it  was  adopted  in  this  constitutive 
modeling  effort.  This  function  has  two  material  parameters:  slope  M  of  the  critical  line, 
and  preconsolidation  pressure  po.  Soils  become  stronger  upon  freezing,  which  is 
characterized  in  the  developed  model  by  increasing  preconsolidation  stress  pa.  At  the 
same  time,  inclinations  of  NCL  and  URL  lines  become  flatter  in  the  u,  In p  space.  A 
reasonable  relative  position  of  normal  compression  lines  for  a  soil  in  both  unfrozen  state 
and  frozen  state  is  shown  in  Figure  21(a).  The  relationship  between  preconsolidation 
stress  po  and  ice  ratio  is  illustrated  in  Figure  21(b).  The  yield  surface  and  its  evolution 
with  ice  ratio  e,  is  illustrated  in  Figure  21(b)  (p0  is  the  preconsolidation  stress  for 
unfrozen  soil,  and  pi  is  the  preconsolidation  stress  for  the  same  soil  in  frozen  state  with 

the  ice  ratio  of  ei  =  eiC ,  /?,'  is  a  reference  stress). 

Consider  freezing  and  loading  path  B-C-D  as  illustrated  in  Figure  21.  A  virgin  unfrozen 
(saturated)  soil  at  state  characterized  by  point  B  is  subjected  to  freezing  under  a  constant 
isotropic  stress.  The  specific  volume  increases  due  to  the  volumetric  expansion  upon 
phase  change.  Then,  isotropic  load  is  added  at  constant  temperature  until  point  D  is 
reached  on  the  isotropic  stress  yield  line  for  the  frozen  soil.  During  the  isotropic  loading 
the  ice  ratio  remains  constant,  ei  =  eiC .  This  is  because  the  ice  ratio  is  related  to  the 

unfrozen  water  content  in  non-segregation  freezing  process,  thus  having  a  unique 
relationship  to  the  temperature. 
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Figure  21.  (a)  Freeze-thaw  thermal  loading  path  in  v-\n(p)  plane;  (b)  “pseudo 
preconsolidation  pressure”  changing  along  BD  due  to  increasing  pore  ice  ratio 

Points  B  and  D  belong  to  the  same  yield  curve  on  p,et  -plane.  Yield  stresses  p({  for 
frozen  states  with  different  ice  ratio  ei  during  freezing  from  point  B  to  C  are  located  on 
this  curve. 

The  specific  volume  at  the  final  point  D  can  be  described  as 

vD  —  vB  +  A  vBC  —  A  vCD  (11) 

where  A vBc  and  A vcd  are  the  increments  due  to  initial  freezing  and  due  to  subsequent 
loading. 

The  expansion  upon  freezing  is  related  to  phase  change  and  is  treated  as  reversible  (it 
reversed  during  thawing).  This  volume  change  can  then  be  calculated  as 

Av/=0.09  et  (12) 

Substituting  equations  (9),  (10),  and  (12)  into  (11),  one  obtains 
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(13) 


vi 


A,  ln-^y  =  (v0  -Aln-^)  +  0.09  e  -at.  In  — 


A> 


Po 


Po 


Po 


Po 


where  /?,'  is  the  reference  pressure  to  locate  the  virgin  compression  lines  for  both  soil 

states.  Consequently,  the  increase  in  isotropic  yield  stress  caused  by  the  ice  ratio  increase 
is  found  as 
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Po 

Po 


=  e 


A-fCf 


\ 

Po_ 

\Po) 


(14) 


The  slopes  of  the  NCL  and  URL  of  frozen  soil  are  both  functions  of  ei .  Because  the 
frozen  soil  has  limited  tendency  for  compression  (pores  are  filled  with  ice),  slope  Af 

drops  dramatically  compared  to  that  for  the  unfrozen  soil.  The  following  relationship  is 
postulated 


Af  =Aexp(-aei)  (15) 

where  «is  a  soil  constant,  k f  is  assumed  to  be  constant  during  freezing,  i.e.,  Kf  =k  . 

In  order  for  equation  (14)  to  yield  po  for  unfrozen  soil  when  e,  =  0,  we  postulate  the 
following  linear  law  defining  the  shift  of  vo  to  v/  as  a  function  of  <?, 

Uf  =v0-(/3-O.O9  )e{  (16) 

where  is  a  soil  constant  (to  be  determined  from  calibration).  Equation  (14)  can, 
therefore,  be  written  as 
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Po 


-fie, 


X~Kf 


—  e 


Po 

v  Po  J 


(17) 


The  plastic  volumetric  strain  due  to  mechanical  load  and  the  thermal  process  is  then 
described  in  the  following  expression 

a8) 

v  Po  v  Po 

The  yield  function  in  the  model  is  similar  to  that  in  Cam  clay  model,  but  with  the 
preconsolidation  pressure  po  dependent  on  ice  ratio  e. 


221  P  a  g  e 


f  =  q2-M2p[p0(ei)-p\  =  0 


(19) 


An  associate  flow  rule  is  assumed  to  describe  plastic  deformation  (i.e.,  the  plastic 
potential  g  =  0  had  the  same  form  as/=  0. 

Thawing  of  a  frozen  soil  decreases  the  amount  of  ice  in  pores,  reduces  “ice  cementation” 
and,  consequently,  reduces  the  strength.  For  non-segregation  frozen  soils,  the  strength 
after  complete  thawing  tends  to  move  back  to  the  same  envelope  as  before  freezing. 
There  will  be  no  further  weakening  since  no  thermally  induced  mass  migration  occurred 
during  freezing,  and  the  change  of  the  soil  fabric  during  the  thermal  process  was  not 
significant. 

5.2  Thermal-mechanical  load  process 

To  illustrate  the  model,  a  freeze-thaw  cycle  along  with  an  external  load  applied  is  shown 
in  Figure  22.  The  thermal  process  is  depicted  with  the  blue  (freezing)  and  red  (thawing) 
lines,  whereas  the  mechanical  loading-unloading  is  depicted  by  the  green  lines. 

A  saturated  soil  specimen  is  preconsolidated  under  isotropic  compression  from  point  A  to 
B.  Then,  under  constant  stress,  the  freezing  process  takes  place  (B-C),  and  the  ice  ratio  at 
point  C  is  eiC .  In  this  process,  the  isotropic  yield  stress  (apparent  preconsolidation  stress 
for  corresponding  frozen  soil)  increases  from  Bi  to  Di  following  equation  (14),  and  it  has 
a  value  of  at  point  Di  (the  soil  has  the  same  apparent  preconsolidation  stress  after  it 

had  been  frozen  at  point  Ci).  The  below-freezing  temperature  is  then  maintained,  and 
isotropic  compression  is  increased  to  reach  the  normal  compression  line  for  frozen  soil  at 
pressure  p °  ;  the  load  is  then  continued  along  the  NCL  to  reach  p(f  .  During  this  loading, 

the  void  ratio  e  is  changing  while  the  ice  ratio  remains  constant.  This  is  based  on  an 
assumption  that  ice  and  soil  skeleton  are  both  incompressible.  From  C  to  D,  the  frozen 
soil  experiences  elastic  behavior,  whereas  from  D  to  E  it  behaves  plastically. 
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Figure  22.  Freeze-thaw  cycle  and  load  path  with  corresponding  yield  curves 

Thereafter,  the  soil  is  thawed  at  constant  stress  p .  Once  ice  starts  melting,  the  soil  can 

no  longer  sustain  load  p  =  ,  and  the  process  of  consolidation  will  start,  moving  the  soil 

to  the  normal  compression  line  for  the  unfrozen  soil  (point  F).  Unloading  from  F  results 
in  an  elastic  rebound  along  the  URL  for  unfrozen  soil  to  G. 


5.3  Implementation  and  application  of  the  model 

The  constitutive  model  was  implemented  into  the  finite  element  system  ABAQUS  using 
subroutine  UMAT  and  UEXPAN  to  solve  boundary  value  problems. 

The  model  was  calibrated  based  on  the  tests  done  by  Qi,  et  al.  (2010).  Parameters 
obtained  from  curve  fitting  into  the  test  data  are  listed  in  Table  5.  The  corresponding 
preconsolidation  stress  curve  with  freezing  temperature  is  shown  in  Figure  23. 
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Figure  23.  Calibration  of  preconsolidation  stress  for  frozen  soil  vs  temperatures 
Table  5  Parameters  and  initial  values  in  simulation 
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The  parameters  A  and  k  are  the  slopes  for  NCL  and  URL  lines  in  compression  plane.  M  is 
the  slope  for  the  critical  state  line.  ps0‘ar>  and  eo  are  the  initial  values  for  the 

preconsolidation  pressure  and  the  void  ratio,  w*,  w ,  and  a  are  parameters  for  unfrozen 
water  content  function. 

To  illustrate  the  response  of  the  model  to  loading  and  thermal  changes,  a  soil  column 
subjected  to  both  the  mechanical  load  and  a  thermal  process  was  simulated.  The  soil 
column  is  0.07  m  in  height.  The  vertical  walls  of  the  column  are  adiabatic  and  rigid.  The 
initial  uniform  temperature  is  1°C,  and  the  initial  vertical  and  horizontal  compressive 
stresses  are  20  kPa  and  10  kPa,  respectively.  Initial  and  boundary  conditions  in  terms  of 
load  and  temperature  at  the  top  of  the  column  are  shown  in  Figure  24.  The  bottom  of  the 
column  is  fixed  and  the  temperature  is  maintained  at  1°C  throughout  the  process.  The 
temperature  distribution  at  t  =  0  and  t  =  300  h  are  shown  in  Figure  25.  More  details  on 
this  simulation  can  be  found  in  Zhang  and  Michalowski  (2013). 
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Figure  24.  Boundary  condition  on  top  of  the  soil  column 

The  relationship  between  void  ratio  and  the  vertical  stress  is  shown  in  Figure  26  for 
element  #18  and  #4.  For  element  #18,  the  temperature  is  about  -4.5°C  when  the  column 
reaches  steady  state  after  freezing,  having  more  pore  ice  and  being  stronger  than  element 
#4  whose  steady  state  temperature  is  around  -0.3°C.  There  is  a  substantial  difference  in 
the  behavior  of  the  two  elements  during  the  load  segment  from  60  to  210  kPa:  while 
element  #18  behaves  elastically,  element  #4  is  elastic  only  until  the  load  reaches  160  kPa, 
and  becomes  elasto-plastic  afterward.  Both  elements  experience  additional  settlement  due 
to  thawing. 
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Figure  25.  Temperature  distribution  at  t  =  0  and  t  =  300  h 


261  P  a  g  e 


Figure  26.  Compression  for  element  #18  (a)  and  element  #4  (b) 


Other  boundary  value  problems  were  also  simulated  using  the  model  developed.  These 
are  not  presented  here,  but  will  be  a  subject  of  future  publications. 

6  Summary  of  the  most  important  results 

The  principal  contribution  of  this  research  is  the  development  of  the  model  for  simulating 
freezing-thawing  cycles  in  frost-susceptible  soils.  The  model  is  capable  of  simulating 
frost  heave  very  effectively,  and  the  authors  view  it  as  the  most  efficient  of  existing 
models.  This  model  also  addresses  the  problem  of  soil  thawing.  This  work  is  a  pioneering 
effort  to  account  for  the  changes  in  soil  strength  associated  with  thawing,  and  predictions 
of  thaw  weakening. 
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